

#  $$$$$$\  $$$$$$$$\       $$$$$$$\             $$\     $$\           
# $$  __$$\ $$  _____|      $$  __$$\            $$ |    \__|          
# $$ /  \__|$$ |            $$ |  $$ | $$$$$$\ $$$$$$\   $$\  $$$$$$\  
# \$$$$$$\  $$$$$\          $$$$$$$  | \____$$\\_$$  _|  $$ |$$  __$$\ 
#  \____$$\ $$  __|         $$  __$$<  $$$$$$$ | $$ |    $$ |$$ /  $$ |
# $$\   $$ |$$ |            $$ |  $$ |$$  __$$ | $$ |$$\ $$ |$$ |  $$ |
# \$$$$$$  |$$$$$$$$\       $$ |  $$ |\$$$$$$$ | \$$$$  |$$ |\$$$$$$  |
#  \______/ \________|      \__|  \__| \_______|  \____/ \__| \______/ 
                                                                     
                                                                    
                                                                     
load("metadata/iv_replicate.rds")
library("latex2exp")
se.ratio <- d$iv_analy_se / d$ols_analy_se
rho <- d$rho
d$name[which(se.ratio < 0.9)]

x0 <- c(0.5, 1, 3, 10, 30, 70)
pdf("graphs/Figure3a.pdf")
xlab <- TeX(r'($\hat{SE}(\hat{\tau}_{2SLS}) / \hat{SE}(\hat{\tau}_{OLS})$)')
par(mar = c(4, 4, 1, 1))
hist(log10(se.ratio),
  breaks = 30, main = "", ylim = c(0, 10), xlim = c(-0.5, 1.9),
  xlab = xlab, ylab = "#Papers", cex.lab = 1.2,
  axes = FALSE)
box()
axis(2)
axis(1, at = log10(x0), labels = x0)
ratio.mean <- mean(se.ratio)
ratio.median <- median(se.ratio)
abline(v = log10(ratio.mean), lty = 2, col = 4, lwd = 2)
abline(v = log10(ratio.median), lty = 2, col = 2, lwd = 2)
text(1.6, 9.8, paste("   Mean =", round(ratio.mean, 1)), cex = 1.3)
text(1.6, 9.3, paste("Median =", round(ratio.median, 1)), cex = 1.3)
graphics.off()


pdf("graphs/Figure3b.pdf", width = 7, height = 7)
xlab <- TeX(r'($|\hat{\rho}(d, \hat{d})|$)')
ylab <- TeX(r'($\hat{SE}(\hat{\tau}_{2SLS}) / \hat{SE}(\hat{\tau}_{OLS})$)')
col <- rep("#77777750", nrow(d))
par(mar = c(4, 4, 1, 1))
plot(1,
  cex = 0.8, type = "n",
  xlim = c(0, 1), ylim = c(-0.4, 1.9),
  xlab = "", ylab = "", axes = FALSE)
abline(h = 0, col = "gray", lty = 2, lwd = 2)
points(rho, log10(se.ratio), col = col, pch = 16, cex = 1.2, lwd = 0.8)
points(rho[45], log10(se.ratio[45]), col = "gray40", pch = 16, cex = 1.2, lwd = 0.8)
points(rho, log10(se.ratio), cex = 1, lwd = 0.8, pch = 1)
box()
y0 <- c(0.5, 1, 3, 10, 30, 70)
axis(1); axis(2, at = log10(y0), labels = y0)
mtext(ylab, 2, 2.5)
mtext(xlab, 1, 2.5)
graphics.off()





# $$$$$$$\        $$\    $$\          $$\                               
# $$  __$$\       $$ |   $$ |         $$ |                              
# $$ |  $$ |      $$ |   $$ |$$$$$$\  $$ |$$\   $$\  $$$$$$\   $$$$$$$\ 
# $$$$$$$  |      \$$\  $$  |\____$$\ $$ |$$ |  $$ |$$  __$$\ $$  _____|
# $$  ____/        \$$\$$  / $$$$$$$ |$$ |$$ |  $$ |$$$$$$$$ |\$$$$$$\  
# $$ |              \$$$  / $$  __$$ |$$ |$$ |  $$ |$$   ____| \____$$\ 
# $$ |               \$  /  \$$$$$$$ |$$ |\$$$$$$  |\$$$$$$$\ $$$$$$$  |
# \__|                \_/    \_______|\__| \______/  \_______|\_______/ 
                                                                     


load("metadata/iv_replicate.rds")
names(d)
sellnull <- which(d$selling_null == 1)
weakiv <- which(d$f_effective < 10)

d$cluster <- ifelse(is.na(d$f_cluster) == TRUE, 0, 1)
(cl <- which(d$cluster == 1))

# p value
d$iv_report_z <- d$iv_report_coef / d$iv_report_se
d$iv_report_p <- 2 * (1 - pnorm(abs(d$iv_report_z)))
d$name[which(d$iv_report_p > 0.05)]

## plot
pdf("graphs/Figure4a.pdf", height = 6, width = 6)
par(mar = c(4, 4, 1, 1))
plot(1,
  type = "n", xlim = c(0, 1), ylim = c(0, 1), axes = FALSE,
  xlab = "Original P-Value", ylab = "P-Value (Bootstrap-c)"
)
box()
abline(v = sqrt(0.05), h = sqrt(0.05), col = 4, lty = 2)
abline(v = sqrt(0.1), h = sqrt(0.1), col = "gray50", lty = 2)
abline(a = 0, b = 1, col = 2, lty = 2)
points(sqrt(d$iv_report_p[-sellnull]), sqrt(d$iv_bootc_p[-sellnull]), pch = 17, col = "#333333")
points(sqrt(d$iv_report_p[sellnull]), sqrt(d$iv_bootc_p[sellnull]), pch = 16, col = "#333333", cex = 1.2)
axis(1, at = sqrt(x0), labels = x0)
axis(2, at = sqrt(x0), labels = x0)
graphics.off()



pdf("graphs/Figure4b.pdf", height = 6, width = 6)
par(mar = c(4, 4, 1, 1))
plot(1,
  type = "n", xlim = c(0, 1), ylim = c(0, 1), axes = FALSE,
  xlab = "Original P-Value", ylab = "P-Value (Bootstrap-t)"
)
box()
abline(v = sqrt(0.05), h = sqrt(0.05), col = 4, lty = 2)
abline(v = sqrt(0.1), h = sqrt(0.1), col = "gray", lty = 2)
abline(a = 0, b = 1, col = 2, lty = 2)
points(sqrt(d$iv_report_p[-sellnull]), sqrt(d$iv_boott_p[-sellnull]), pch = 17, col = "#333333")
points(sqrt(d$iv_report_p[sellnull]), sqrt(d$iv_boott_p[sellnull]), pch = 16, col = "#333333", cex = 1.2)
# points(sqrt(d$iv_report_p[weakiv]), sqrt(d$iv_boott_p[weakiv]), pch = 17, col = 2)
axis(1, at = sqrt(x0), labels = x0)
axis(2, at = sqrt(x0), labels = x0)
graphics.off()



pdf("graphs/Figure4c.pdf", height = 6, width = 6)
bounded <- which(d$AR_bounded == TRUE)
par(mar = c(4, 4, 1, 1))
plot(1,
  type = "n", xlim = c(0, 1), ylim = c(0, 1), axes = FALSE,
  xlab = "Original P-Value", ylab = "P-Value (Anderson-Rubin)"
)
box()
abline(v = sqrt(0.05), h = sqrt(0.05), col = 4, lty = 2)
abline(v = sqrt(0.1), h = sqrt(0.1), col = "gray50", lty = 2)
abline(a = 0, b = 1, col = 2, lty = 2)
points(sqrt(d$iv_report_p[-sellnull]), sqrt(d$AR_p[-sellnull]), pch = 2, cex = 0.9, col = "#333333")
points(sqrt(d$iv_report_p[setdiff(bounded, sellnull)]), sqrt(d$AR_p[setdiff(bounded, sellnull)]), pch = 17, col = "#333333")
points(sqrt(d$iv_report_p[sellnull]), sqrt(d$AR_p[sellnull]), pch = 16, col = "#333333", cex = 1.2)
# points(sqrt(d$iv_report_p[weakiv]), sqrt(d$AR_p[weakiv]), pch = 17, col = 2)
axis(1, at = sqrt(x0), labels = x0)
axis(2, at = sqrt(x0), labels = x0)
graphics.off()


#   $$\      $$$$$$\                $$$$$$$$\ 
#   $$ |    $$  __$$\               $$  _____|
# $$$$$$\   $$ /  \__|     $$$$$$$\ $$ |      
# \_$$  _|  $$$$\ $$$$$$\ $$  _____|$$$$$\    
#   $$ |    $$  _|\______|$$ /      $$  __|   
#   $$ |$$\ $$ |          $$ |      $$ |      
#   \$$$$  |$$ |          \$$$$$$$\ $$ |      
#    \____/ \__|           \_______|\__|      


library(ivDiag)
library("latex2exp")

F0 <- c(seq(2, 10.3, 0.1), 100000)^2
cF0 <- c(
  18.66, 9.74, 7.37, 6.18, 5.43, 4.92, 4.54, 4.25, 4.01, 3.82, 3.65, 3.51, 3.39, 3.29, 3.19, 3.11, 3.03,
  2.97, 2.91, 2.85, 2.80, 2.75, 2.71, 2.67, 2.63, 2.60, 2.57, 2.54, 2.51, 2.48, 2.46, 2.43, 2.41, 2.39, 2.37,
  2.35, 2.33, 2.32, 2.30, 2.29, 2.27, 2.26, 2.24, 2.23, 2.22, 2.21, 2.20, 2.19, 2.17, 2.16, 2.16, 2.15, 2.14,
  2.13, 2.12, 2.11, 2.10, 2.10, 2.09, 2.08, 2.08, 2.07, 2.06, 2.06, 2.05, 2.04, 2.04, 2.03, 2.03, 2.02, 2.02,
  2.01, 2.01, 2.00, 2.00, 1.99, 1.99, 1.99, 1.98, 1.98, 1.97, 1.97, 1.97, 1.96, 1.96
)
F0s <- F0 / (10 + F0)
cF0s <- cF0^2 / (1 + cF0^2 / 1.96^2)


load("metadata/iv_replicate.rds")
d <- d[which(d$p_iv == 1), ]
dim(d) # 59
tstat <- abs(d$iv_coef / d$iv_analy_se)
round(tstat, 3)
good <- which(tstat >= d$tF_cF)
bad <- which(tstat < d$tF_cF)
length(good) # 42
length(bad) # 17
length(bad) / nrow(d) # 29%
nrow(d)
d$name[which(tstat < d$tF_cF & tstat > 1.96)]

pdf("graphs/Figure4d.pdf", height = 6, width = 6)
par(mar = c(4, 4, 1, 1))
x <- d$f_effective / (10 + d$f_effective)
y <- tstat^2 / (1 + tstat^2 / 1.96^2)
plot(x, y,
  xlim = c(0.2, 1.01), ylim = c(0, 3.75),
  axes = FALSE, xlab = "", ylab = "", type = "n"
)
box()
mtext("First Stage Partial F Statistic (rescaled)", 1, 2.5)
mtext(TeX(r'($t^2$ Statistic for 2SLS Estimates (rescaled))'), 2, 2.5)
lines(F0s, cF0s, col = 4, lwd = 2)
points(x[good], y[good], col = "green", pch = 17)
points(x[setdiff(bad, sellnull)], y[setdiff(bad, sellnull)], col = "pink", pch = 17)
points(x[-sellnull], y[-sellnull], pch = 2)
points(x[sellnull], y[sellnull], col = "pink", pch = 16, cex = 1.2)
points(x[sellnull], y[sellnull], pch = 1, cex = 1.2)
v0 <- c(3.84, 10, 104.7)
abline(v = v0 / (10 + v0), col = "gray50", lty = 2)
axis(1, at = v0 / (10 + v0), labels = expression(1.96^2, 10, 104.7))
h0 <- c(1, 1.96, 3.43)
abline(h = h0^2 / (1 + h0^2 / 1.96^2), col = "gray50", lty = 2)
axis(2, at = h0^2 / (1 + h0^2 / 1.96^2), labels = expression(1^2, 1.96^2, 3.43^2))
text(0.53, 3, "c(F)", col = 4, cex = 1.2)
graphics.off()





